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Hemorrhoidal artery ligation (HAL) has become universally accepted 
minimally invasive treatment of hemorrhoids disease. HAL involves precise 
identification of the superior rectal arteries supplying hemorrhoidal tissues 
using ultrasonic Doppler principles. During this process, at least there are 
three distinct sounds may be encountered by the surgeon. Only the pulsing 
Doppler sound is useful as it indicates the presence of hemorrhoidal artery. 
The accuracy based on traditional auscultation is commonly affected by 
surgeon’s hearing sensitivity and clinical experience. Therefore, automatic 
Doppler blood flow sound will be a great help in locating hemorrhoidal 
arteries. In this paper, a method based on the center frequency and kurtosis 
features extracted from Burg’s power spectral density (PSD) to distinguish 
three different types of Doppler blood flow sound signal during HAL 
procedure is proposed. Separability measurement was carried out using K- 
means clustering with the city block distance and three clusters 


corresponding to different sound types are successfully formed. In terms of 
arterial sound detection, an accuracy of 94.11% can be achieved. This result 
suggests that centre frequency, kurtosis, and maybe some other statistical 
features extracted from Burg PSD have the potential to be utilized as a 
means in automatic Doppler blood flow sound recognition. 
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1. INTRODUCTION 

The surgical intervention techniques of hemorrhoids disease have advanced over past the three 
decades. Traditional surgical hemorrhoidal excision is notorious for significant postoperative pain and 
impaired anal function due to the deformation of anal canal anatomy. These problems have encouraged 
surgeons to adopt minimally invasive procedures in treating hemorrhoids disease. Since its invention in 1995 
by Morinaga et al. [1], hemorrhoidal artery ligation (HAL) has become a widely accepted procedure for 
handling hemorrhoids disease. In comparison with excisional methods, HAL excels at the absence of anal 
lacerations, which in turn reduces patient’s postoperative pain and uncomfortable feelings. Anal anatomy and 
physiology preservation cannot be underrated as well [2]. Several studies suggest that early and midterm 
outcomes have shown great recovery and positive patient feedbacks [3]. Dearterialization technique, in 
comparison with other non—excisional procedures, such as stapled hemorrhoidopexy may have additional 
benefit of minimizing life-threatening complications [4]. 

The reasoning for HAL procedure is based on the presumption that blood flowing through rectal 
artery to hemorrhoidal tissues cannot recirculate given the absence of capillary intercession between the 
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arterial and venous systems inside the anal canal. This blood overflow eventually causes the dilatation of 
hemorrhoidal tissues. HAL involves the precise marking of of rectal artery branches supplying hemorrhoidal 
tissues using an ultrasonic transducer linked to the outer surface of a proctoscope. Identified hemorrhoidal 
arteries, which are commonly surrounding the anus in the odd-numbered clock positions [5], are then suture 
ligated through a ligation window within the proctoscope. Ligation of these arteries significantly reduces 
blood passage into the hemorrhoidal venous plexi and later causes the hemorrhoidal tissues shrink gradually. 

Hemorrhoidal arteries identification in HAL procedure is based on the Doppler ultrasound 
principles. To properly detect each hemorrhoidal artery within the rectal column, the surgeon conducting the 
HAL surgery relies solely on the Doppler sound emitted by the Doppler processor. There are two major 
constrains of Doppler sound auscultation. The first one is related to the interpretation of the audible output by 
the surgeon, as pointed out by Mowatt et al [6]. The inaccurate interpretation of Doppler ultrasound outputs 
may arise from inadequate training, lack of experience, and time constraints [7]-[9]. The second one is that 
the surgeon must listen very carefully to detect subtle changes in Doppler ultrasound outputs [10]. 

Generally speaking, the accuracy based on traditional auscultation is highly affected by observer’s 
hearing sensitivity and clinical experience [11]. Therefore, a computerized assisting tool to analyze Doppler 
ultrasound signal can be beneficial to recognize certain Doppler sound that indicates haemorrhoidal arterial 
blood flow. However, Doppler sound signal produced by ultrasound machine is rarely recorded and 
investigated, including the one associated with hemorrhoidal artery [12]. Thus far, most research has focused 
on other biological sounds such as heart sound (HS) and lung sound (LS). Therefore adopting methods from 
studies regarding HS and LS analysis is very relevant to keep methods that will be used right on the cutting 
edge. Doppler blood flow sound, HS, and LS share common audible frequency spectrum and they all 
represent a non-stationary signal [10], [13]-[17]. 

Despite the limited number of studies addressing Doppler sound signal analysis, some literatures are 
accessible that could serve as a steppingstone. One of the preliminary works in investigating Doppler sound 
signal characteristics was conducted by Guler and Kara [18] who used fast fourier transform (FFT) and 
autoregressive (AR) modelling methods to perform spectral analysis of Doppler sound signal produced by 
healthy and stenosis mitral valve. It was reported that spectrogram of healthy subject had two peaks and a 
valley between them while for the subject suffering mitral stenosis, the valley was absent. Later, Guler et al. 
Proposed an automated diagnosis method to distinguish Doppler sound signals acquired from both healthy 
subjects and subjects having ophthalmic arterial (OA) and internal carotid arterial (ICA) diseases [19]. The 
method involved statistical features: maximum, mean, minimum, and standard deviation calculated from time 
— frequency representations of Doppler blood flow sound signals. The extracted features of different 
categories of OA and ICA signals were different from each other so they can serve as discriminating 
parameters in classifying Doppler signals. In this study, the authors also mentioned that performing spectral 
analysis is the most complete way to present blood flow information contained in the Doppler shift signal. 
Maximum envelope of the carotid artery Doppler spectrogram derived from AR method was used by Ozsen 
for diagnosing atherosclerosis disease [20]. It was found that subjects with atherosclerosis have lower 
maximum frequency envelope values. According to Tedim ef al. [21] was the first to search automatic 
methods to improve venous air embolism (VAE) detection using precordial Doppler placed on human 
patients during real surgeries. VAE may occur during surgical procedures and may lead to serious 
complications such as stroke or cardiovascular failure. VAE event, which is indicated by turbulent blood 
flow, was simulated by using saline administrations through peripheral and central catheters. Frequency 
component of Doppler sound signal was examined by power spectral density (PSD) estimation by using 
Welch method. Different kinds of spectral features including frequency correlating with maximum power of 
PSD, 95% power of PSD, and half power of PSD were then extracted to characterize the power distribution 
attribute. Observations showed that frequency corresponding to 95% power of PSD was the most sensitive to 
detect blood flow turbulence whereas frequency corresponding to half power of PSD remained unchanged. 

As for HS and LS, more research had been conducted to study their characteristics [22]. A research 
performed by Theo et al. [23]. Assessed the feasibility of using the analysis of power spectral to compute the 
ratio between the sum of the powers in the N highest peaks in the power spectral to that of the entire power 
spectralfor each fundamental HS state. This feature can be used to classify normal or abnormal HS. Kristomo 
et al. [24] presented a method based on the statistics calculation from autoregressive—power spectral density 
(AR-PSD) extracted from HS. These statistics were used as inputs for the classification of abnormal HS into 
9 types. Research by Schmidt et al. [25] weak murmurs in HS were used to detect coronary artery disease 
(CAD) by dividing HS into five different frequency bands and nine features were extracted from those bands. 
A multivariate classifier based on quadratic discriminant analysis (QDA) was used to construct CAD score 
from the features. According to Akanksha et al. [26] with a similar goal collected HS from four locations on 
chest and was analysed using cross power spectral density (CPSD). Four separate features, including relative 
power, power ratio to adjacent subbands, mean and standard deviation were computed from CPSD. Support 
vector machine (SVM) was then used as a classifier. Different power spectral characteristics estimated from 
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LS may also be used to indicate pulmonary chronic lung diseases, as investigated in [27], [28]. K—means 
algorithm, one of the most widely used clustering approaches, was also used to separate Mel-frequency 
spectral coefficient (MFCC) features within first and second HS segments (S1 and S2, respectively) into two 
groups [29]. A population center vector for each group was then calculated to form a supervector, which was 
used as input for a deep neural network (DNN) classifier. The results revealed improvement in terms of 
accuracy by employing k-means prior to the classification process. Jamal et al. [30] added murmur 
components in HS analysis in addition to S1 and S2 detection by using time—frequency method. Research by 
Hadiyoso et al. [31] took advanced approach by implementing an electronic stethoscope that enabled the 
findings on auscultation to be sent to a cloud server for more advanced analysis. Android—based mobile 
phone and Steder application were used as a communication device and sound data manager, respectively. 
The prototype of the system was successfully developed and usable but had not been tested yet in real 
situations. 

From the literature review, it is discovered that power spectral analysis has promising opportunity in 
characterizing Doppler sounds. In this paper, a simple method to distinguish different types of Doppler blood flow 
sound signal during HAL procedure is proposed. Burg’s PSD is applied to estimate the band power of different 
blood flow sound signals which is grouped into three classes. Statistical features extracted from the signal power 
bands are expected to provide means to accurately distinguish the sound signals into predetermined classes. 
Performance evaluation is based on t-test and k-means clustering with city block distance. 


2. MATERIALS AND METHODS 
2.1. Processing steps 

The standard procedure for computerized PCG processing and analysis can be outlined into the 
following sequences: i) pre-processing; ii) feature extraction; iii) classifier modeling [16]. The characterization 
problem in this study was also resolved according to the procedure with several required adjusments. 

Doppler sounds were acquired using an 8 MHz transducer that gently moved around the wrist area 
to scan for any blood movement inside radial artery. Rather than a real site inside the anal canal, the wrist 
area was selected because of the ease of access and similar Doppler produced. This approach was approved 
by collaborating surgeon. The processor received the frequency-shifted signal from the transducer and 
generated a Doppler signal with a frequency within the audible band. A Bluetooth-enabled portable amplifier 
was used to amplify the signal. The Doppler audible sounds were captured by using a digital sound recorder 
with a 44.1 kHz sampling rate. 

The first processing step was to reduce the sound sample to a single arterial pulse cycle, with 
approximately 0.7 s duration. Then the sample was resampled at a sampling frequency of 8820 Hz and 
followed with amplitude normalization. Next, the power spectral of the sample was analyzed using Burg 
method while transforming the sample into a frequency domain. After that, features were extracted from PSD 
curve and analyzed for the most suitable features for classification. Finally, separability measurement was 
performed to determine whether selected features can be used to discriminate one type of sound from 
another. A total of 51 sound samples used, separated into classes: arterial, venous, and probe rubbing sounds. 
Figure 1 shows the system flowchart for the characterization of Doppler blood flow sound during HAL. 
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Figure 1. Flowchart of characterization process 


2.2. Doppler blood flow detection principle 

The Doppler effect refers to the frequency shift of a signal reflected from a moving object to its 
observer perception. Doppler signal can be possibly obtained from blood flow since blood contains red blood 
cells as scattering particles. If the blood flow is toward the transducer, the perceived frequency will be higher 
and vice versa. Figure 2 illustrates an ultrasound beam being transmitted from the transducer toward red 
blood cells at a certain angle of attack 6. The beam hits the objects and returns with a Doppler-shifted 
frequency to the receiver. The Doppler shift fd of an ultrasound signal with the nominal frequency fc is given 


by (1): 
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2v cos 0 


p= (1) 


Where v is the velocity of the scattering objects and c is the propagation speed of sound waves in the 
medium. The scattering object in this case is red blood cells flowing inside haemorrhoidal vessels with a 
normal velocity of 20 to 750 mm/s [32]. Sound waves travel in soft tissue at an average speed of 1540 m/s 
[33]. The frequency f, for blood flow detection application is in the range of 2 to 10 MHz, resulting in a 
Doppler shift in the audible range of 0.2 to 7.5 kHz. 


blood vessel 
Oad OO 


Figure 2. Doppler—shifted frequency fd appears as an audible sound as the reflected beam received by the 
transducer 


Doppler sound waveforms are commonly defined as triphasic. The key components of a Doppler sound 
waveform include: (i) primary forward flow due to contraction of left ventricle during systole; (ii) transient reversal 
flow as a result of backscattering from outflow bed with a high resistance in early diastole; (iii) secondary forward 
flow due to reflection from a closed aortic valve during late systole [34]. In Figure 3, both Doppler sound and PCG 
waveform are plotted together to show the relationship between them in time domain. 
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Figure 3. Doppler sound waveform and four states of cardiac cycle (S1, systole, S2, diastole) presented [34], [35] 


Different Doppler signals are produced depending on the position of the artery, distance from the 
transducer, and direction of blood flow in relation to the ultrasound waves emitted by the transducer. In terms 
of intensity, the more perpendicular the blood flow to the ultrasound waves the higher the Doppler signal. 
Conversely, the more parallel the blood flow is, the lower the signal. The hemorrhoidal arteries are sought by 
slowly rotating and tilting the proctoscope inside the anal canal. During this process, at least three different 
sounds may be encountered by the surgeon, namely; pulsing, wheezing, and thumping sounds. Pulsing sound 
reflects pulsatile blood movement inside an artery, wheezing sound reflects non-pulsatile venous velocity, 
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and thumping sound results from rubbing action between transducer and rectum wall. Only the pulsing 
Doppler sound indicates the presence of hemorrhoidal artery, characterized by middle—-low frequency 
spectrum (up to 4 kHz) with pulsing amplitude. The return echoes from steady blood stream may occur when 
venous blood vessels pass through the region of Doppler sensitivity [36]. 


2.3. PSD estimation 

PSD describes how the power of a signal is distributed along frequency. The Wiener—Khintchine 
theory suggests that the PSD of a wide—sense stationary random process x(t) is the Fourier transformation of 
the autocorrelation of the signal [37]. This definition can be derived using (2) to (4): 


F[x(t)] = X(@) = [2 x(@e/'dt (2) 
R(t) = x(t) * x(=T) = J" x(e)x(t + 1) dt (3) 
F[R(t)] = fo, Re“ dt (4) 


Autoregressive (AR) based model is one of the parametric methods used to estimate the PSD of a 
signal. AR models use a linear predictive method to extrapolate the signal beyond its existing values and 
predict the output of a system based on previous outputs [37]. Such parametric method provides a fine 
frequency resolution but is sensitive to order selection. A model with an extremely low order produces an 
over—smoothed spectral whereas an extremely order tends to produce false low-level peaks in the spectral. 
An AR model can be defined as in (5): 


x(t) = ro aiXti + E(t) (5) 


Where a; are the AR coefficients, p is the order of model, and e(t) represents noise. The PSD estimation by 
using an AR model is expressed by (6): 


a E- 

P, = —— 2 6 

ROS a © 
Burg is one of the several methods available to estimate AR coefficients. In Burg method, AR parameters are 
found by minimizing (least squares) errors of sums of forward and backward linear prediction. Such 
minimization occurs with the AR coefficients constrained to satisfy the Levinson—Durbin recursion 
algorithm. 


2.4. Statistical features 
Abundant statistical features are available to be extracted from data both in time and frequency 
domains. Some of the features that are relevant to this study are listed below. 
a. Mean: the arithmetic mean p is the average of values {x1 X2,..., Xm} located within a certain data range. 
It is expressed by (7): 


— tym 


B= 7 Lisi Xi (7) 


b. Standard deviation: standard deviation can be calculated using (8). It is an indication how values 
{X1, X2...» Xm} are spread out. 


o= [ASR —w? (8) 


c. Skewness: skewness is one of the features that is used to measure the shape of the data set, particularly 
the asymmetry. It is expressed as in (9): 


Sk = £ (9) 


o* 


Where u; is the 3‘ moment about the mean. 
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d. Kurtosis: similar to skewness, kurtosis provides information about the shape of the data set. It is a 
measure of the “tailedness” of the probability distribution of the data. The formula is given by (10): 


Ku =" (10) 


ot 


e. Entropy: the entropy of a random variable is used to determine how unpredictable it is. The entropy H of a 
discrete variable X with possible values {x,, x2,..., Xm } and probability mass function p(X) is given by (11): 


H(X) = — Wiz, P (x;). log, P (x;) (11) 


f. Frequency related: Figure 4 shows the typical PSD estimation result in the frequency domain and 
graphical parameters that can be extracted from it. 
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Figure 4. Typical PSD estimation curve with available features 


2.5. K—means clustering 

Clustering is one of the most popular analytic data identification methods used to provide a view of 
the the data structure. This process identifies homogenous subgroups within the data such that data points in 
each group are as similar as possible in accordance with several similarity measures. Clustering is 
categorized as unsupervised learning given that no ground truth exists for evaluation of the clustering 
performance. 

K—means clustering is considered as one of the most usable clustering algorithms since it is simple and 
fast. This procedure iteratively attempts to partition a certain dataset into k pre-defined individual non-overlapping 
clusters, in which each datum belongs only to one cluster. It works through the following steps [38]: 

a. Choose k starting centers arbitralily C = {c,,..., cx}. 

b. Set cluster C; to be the set of points in X that are closer to c; than they are to c; for all j + i, for each 
i {1,..., k}. 

c. Set c; to be the center of mass of all points in C: c; = req Ërec: x, for each i E {1,...,k}. 


d. Repeat steps 2 and 3 until C remains constant. 
For the k-means problem, given an integer k and a set of n data points X C R®. It is required to choose k 
centers C so as to minimize the potential function in (12) [39]: 


Q = Yxec MiNcec ||x — cll? (12) 


In k-means clustering, data points are assigned to the centroid with a minimum distance based on 
the values found. Therefore, distance calculation plays an important role in this clustering algorithm. Several 
different techniques are available to compute distance between two points such as city block, Euclidean, 
cosine, and correlation distances. For city block distance between two points a and b with k dimensions is 
defined as in (13) [39]: 


d(a,b) = Y*_,|a;_b;| (13) 
Separability evaluation was performed to observe whether statistical features extracted from PSD 
can be used to differentiate one type of sound from another. K—means clustering with a city block distance 


was used for the testing purpose. Clustering performance can be evaluated using a measure such as accuracy 
which is defined as in (14): 


Burg power spectral density-based characterization of Doppler blood flow sound ... (Daniel Santoso) 


1338 O ISSN: 2302-9285 


truepos-ttTueneg 


accuracy = x100% (14) 


trueposttrueneg+falsepostfalseneg 


3. RESULTS AND DISCUSSION 

Three kinds of sounds are normally perceived by the surgeon from the Doppler processor through 
the HAL procedure. However, only the pulsing blood flow sound is useful for the surgeon as a location 
marker of hemorrhoidal arteries. The wheezing sound is produced when the ultrasound wave hits blood 
particles flowing in the vein. When proctoscope is being rotated inside the anal canal in attempt to search 
hemorrhoidal arteries, random rubbing sound is also produced by Doppler processor. The latter two 
mentioned sounds are not sounds of interest for the surgeon. Figure 5 shows the three kinds of Doppler 
sounds that commonly encountered during HAL procedure in the time domain. 
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Figure 5. Sound fragments representing the three classes of Doppler blood flow sound 


Using Burg’s PSD the Doppler sounds were transformed into the frequency domain after 
normalization and downsampling. The purpose of this transformation is to obtain the features of each 
Doppler sound mentioned previously. The transformation results are illustrated in Figure 6. Figures 6(a), (b), 
and (c) are the results of arterial, venous, and rubbing sound sample respectively. A brief observation of the 
figure reveals that all three Doppler sounds have significant frequency components up to 2000 Hz. However, 
different sounds exhibit varied peak characteristics. The PSD estimation for arterial sound has high-low 
peaks while venous sound counterpart has two nearly identical peaks. In the other hand, rubbing sound PSD 
estimation only consists of a single peak with narrower bandwith. In this study, only features related to shape 
of data set were extracted from PSD estimation that is center frequency, kurtosis, entropy, and skewness. 
Table 1 shows the PSD estimation features extraction result. 


Table 1. Extracted features from PSD estimation of all Doppler blood flow sounds 


Doppler sound Center frequency (Hz) Kurtosis Entropy Skewness 
Arterial 434.71 + 17.79 26.37£7.07 2.9840.15 4.48 +0.70 
Venous 486.39 + 31.37 13.0744.17  2.58+0.57 3.13 + 0.58 
Rubbing 314.13 + 14.32 55.35 + 8.31 1.38+0.14 6.91 +0.55 


From this data, it can be observed that the center frequencies of all sounds are relatively uniformed, 
indicated by a slight standard deviation. For kurtosis feature, despite the wide deviation, no overlap is 
observed between sounds. The entropies of all sounds converge into small values with a relatively large 
standard deviation causing the entropies of arterial and venous sounds to overlap. Skewness shows slightly 
higher average values and lower standard deviation compared to entropy but without the overlap. 
Furthermore, t-test was performed to statistically demonstrate class separation, as described in [40]. 
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Figure 6. PSD estimation of (a) arterial, (b) venous, and (c) rubbing sound sample 


Figure 7 illustrates that normalized value statistical features is separable. This result is confirmed by 
t-test tabulated in Table 2. The t-test is an inferential statistic calculated from the value of center frequency, 
kurtosis, entropy, and skewness features to determine whether there is significant difference between the 
means of two groups, arterial or non—arterial (venous and rubbing) sounds in this case. The range of p—value 
was found by calculating t-score with the significance level of 0.05. Table 2 proves that all of the p—values 
fall in the significance range, therefore demonstrating that values of features are separated. 
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Figure 7. Normalized value of statistical features 


In this study, center frequency and kurtosis were chosen as key features to distinguish arterial, 
venous, and rubbing sounds because of their particular characteristics. Center frequency is relatively uniform 
while kurtosis is quite spaced away between one sound type to another. The center frequency and kurtosis as 
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discerning features were tested using k-means clustering with the city block distance. At the end of iterative 
process, the data should be distributed into the three different clusters corresponding to three Doppler sounds 
discussed previously. The result of k-means clustering is illustrated as scatter diagram in Figure 8. 

It can be readily seen that the plotting of two features center frequency and kurtosis shows three 
different areas along with their own center points. The system is expected to correctly detect arterial sound, 
therefore accuracy for this detection needs to be investigated. Based on the clustering, the numbers of true 
positive, true negative, false positive, and false negative samples are 15, 33, 2, and 1 sample respectively. 
This yields accuracy of 94.11% for arterial sound detection by using (14). These findings confirm our 
hypothesis that arterial, venous, and rubbing sounds may be recognized by extracting statistical features from 
PSD estimation of the sound signal. In addition, the fact that power spectral analysis is shown to be effective 
in Doppler blood flow sound characterization is consistent with the discoveries which mainly reported in 
[18]-[21], [23]-[29]. 


Table 2. T-score and p—value calculated from four statistical features 


Feature t-score p-value 
Center frequency 10.4981 0.0129 
Kurtosis -9.5670 0.0104 
Entropy 35.6868 <0.0001 
Skewness -6.6680 0.0416 
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Figure 8. Scatter diagram of center frequency and kurtosis 


4. CONCLUSION 

In this work, an effort has been conducted to characterize and classify three different Doppler 
sounds that usually encountered during HAL procedure. The work involved the Burg method to estimate 
PSD of the sound signal and the estimation curves suggested that each sound type has its own specific 
shape.Statistical features were extracted from the PSD estimation and it was discovered by using t-test that 
center frequency, kurtosis, entropy, and skewness normalized value is separable. Therefore, those features are 
suitable to describe PSD estimation curve quantitatively. Separability measurement was further carried out 
using k-means clustering with center frequency and kurtosis as selected features and three clusters 
corresponding to the three different sound types were successfully formed. In terms of arterial sound 
detection, an accuracy of 94.11% can be achieved without further signal processing. These results suggest 
that center frequency, kurtosis, and maybe some other statistical features extracted from Burg PSD have the 
potential to be utilized as a means in automatic Doppler blood flow sound recognition. In the near future with 
the collection of more data, the effectiveness of the proposed method should be investigated. The ultimate 
goal is real-time hardware implementation of the mentioned system that can be used as an automatic tool in 
clinical setup. 
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